La interpretación de los coeficientes en modelos lineales (generalizados) es más sutil de lo que puedes darte cuenta, y tiene consecuencias en cómo probamos hipótesis y reportamos hallazgos. Comenzaremos hablando de las interpretaciones marginales vs. condicionales de los parámetros del modelo.
En este ejemplo, modelamos la altura de las plantas en función de la altitud y la temperatura. Estas variables están correlacionadas de manera negativa: hace más frío a mayor altitud. Comenzamos simulando algunos datos para reflejar esto.
library(mvtnorm)# Specify the sample sizeN <-1000# Specify the correlation between altitude and temperaturerho <--0.4# This line of code creates two correlated variablesX <-rmvnorm(N, mean =c(0, 0), sigma =matrix(c(1, rho, rho, 1), 2, 2))# Extract the first and second columns to vectors named temp and alt and plottemp <- X[, 1]alt <- X[, 2]plot(temp, alt)
Ahora podemos simular algunos datos de altura de las plantas. Aquí decimos que la altura media de las plantas es 2 (cuando todas las demás variables son 0). A medida que la temperatura aumenta en una unidad (manteniendo la altitud constante), la media de la altura aumentará en 1 unidad (beta[2] = 1), y de manera similar, a medida que aumentas la altitud en 1 unidad (manteniendo la temperatura constante), la media de la altura disminuirá en 1 (beta[3] = -1). La altura sigue una distribución normal con esta media y una desviación estándar de 2.
La interpretación de estos coeficientes es que si mantienes todo lo demás en el modelo constante (es decir, la temperatura) y agregas 1 a la altitud, entonces la altura media estimada disminuirá en 1.09. Ten en cuenta que el coeficiente depende de las unidades en las que se mide la altitud. Si la altitud se mide en metros, entonces el coeficiente te dice qué sucede cuando subes 1 metro.
La intersección es el valor predicho cuando todas las demás variables se establecen en 0, lo cual a veces tiene sentido (aquí sería la altura de las plantas a nivel del mar y a 0 grados de temperatura). Otras veces, 0 no es un valor significativo, y si deseas interpretar la intersección, podría tener sentido reescalar tus otras variables para que su media sea 0. Si haces esto, entonces la intersección es el valor predicho cuando todas las demás variables están en su nivel medio.
¿Y si ahora tuviéramos un modelo solo con la temperatura?
lm1 <-lm(height ~ temp)lm1$coefficients
(Intercept) temp
1.999833 1.467180
El coeficiente de temperatura ahora es 1.38, ¿qué está pasando? La altitud es un predictor importante de la altura de las plantas, y parte de la información sobre la altitud está contenida en la temperatura (recuerda que están correlacionadas, por lo que a medida que la altitud aumenta, la temperatura disminuye). El modelo tiene en cuenta esto cambiando el efecto de la temperatura para tener en cuenta la información que contiene sobre la altitud. Observa que el coeficiente de temperatura está incorrecto en aproximadamente 0.4, que es la cantidad de correlación entre las variables.
Nota: Cuando los estadísticos hablan de esto, usan las palabras condicional y marginal. Condicional es el efecto de una variable cuando las demás se mantienen constantes (como en lm_both), mientras que marginal es el efecto global (como en lm1). Nota: Si utilizas el código anterior para simular tus propios conjuntos de datos, obtendrás valores ligeramente diferentes para los coeficientes.
Pruebas de hipótesis
Esta distinción tiene muchas consecuencias tanto para el modelado como para las pruebas de hipótesis. Generemos algunos datos en los que la altitud predice la altura, y la temperatura no tiene (información adicional), y luego probemos la temperatura.
mu <-2-1* altheight <-rnorm(N, mu, sd =2)mod_temp <-lm(height ~ temp)summary(mod_temp)
Call:
lm(formula = height ~ temp)
Residuals:
Min 1Q Median 3Q Max
-6.5209 -1.4534 -0.0109 1.4639 7.7749
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.99228 0.06953 28.655 < 2e-16 ***
temp 0.46470 0.06656 6.981 5.33e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.197 on 998 degrees of freedom
Multiple R-squared: 0.04656, Adjusted R-squared: 0.04561
F-statistic: 48.74 on 1 and 998 DF, p-value: 5.326e-12
anova(mod_temp)
Analysis of Variance Table
Response: height
Df Sum Sq Mean Sq F value Pr(>F)
temp 1 235.2 235.163 48.739 5.326e-12 ***
Residuals 998 4815.2 4.825
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
La salida de este modelo nos está indicando que hay un efecto de la temperatura, aunque técnicamente no lo haya. No nos está dando información falsa si entendemos cómo interpretar los resultados del modelo. Debido a que la temperatura está correlacionada con la altitud, y hay un efecto de la altitud, cuando la altitud no está en el modelo, el modelo nos dice en general que hay un efecto de la temperatura que aumenta la altura en 0.45 (recuerda que la correlación fue 0.4). Si nuestra hipótesis es “¿La altura de las plantas cambia con la temperatura?”, la respuesta es sí, a mayor temperatura, más altas son las plantas.
Pero, ¿qué pasa con la altitud? Sabemos que el efecto de la temperatura que observamos se debe a que está correlacionado con la altitud, la temperatura no predice directamente la altura. Si queremos saber si hay un efecto de la temperatura después de controlar la altitud (manteniendo la altitud constante, es decir, condicional), entonces ajustamos el modelo con la altitud y luego probamos la temperatura.
mod_temp_alt <-lm(height ~ alt + temp)summary(mod_temp_alt)
Call:
lm(formula = height ~ alt + temp)
Residuals:
Min 1Q Median 3Q Max
-5.3731 -1.3011 0.0368 1.3478 5.9044
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.97272 0.06173 31.960 <2e-16 ***
alt -1.07771 0.06563 -16.422 <2e-16 ***
temp 0.02500 0.06487 0.385 0.7
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.95 on 997 degrees of freedom
Multiple R-squared: 0.2496, Adjusted R-squared: 0.2481
F-statistic: 165.8 on 2 and 997 DF, p-value: < 2.2e-16
anova(mod_temp_alt)
Analysis of Variance Table
Response: height
Df Sum Sq Mean Sq F value Pr(>F)
alt 1 1259.8 1259.80 331.4004 <2e-16 ***
temp 1 0.6 0.56 0.1485 0.7
Residuals 997 3790.0 3.80
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El p-valor es aproximadamente 0.95, por lo que no tenemos evidencia de un efecto de la temperatura después de controlar la altitud.
Nota: La distinción entre interpretaciones condicionales y marginales también es válida para modelos lineales generalizados y modelos mixtos.
Covariables categóricas
Cuando tenemos covariables categóricas (por ejemplo, tratamiento), hay varias formas de codificar el modelo, lo que dará diferentes interpretaciones para los coeficientes. Simulemos 120 puntos de datos con 40 en cada uno de los tres niveles de un tratamiento categórico.
N <-120# The effect of treatmenttrt.n <-rep(c(-1, 0, 1), N /3)mu <-2+1* trt.n# Labels for the treatmenttreatment <-factor(rep(c("low", "med", "high"), N /3)) # group labels# Create, Y, a normally distributed response variable and plot against treatmentY <-rnorm(N, mu, sd =2)boxplot(Y ~ treatment)
Si introducimos el tratamiento como una covariable de la forma habitual, el modelo elegirá un tratamiento de referencia (aquí será “high” ya que los niveles se ordenan alfabéticamente), de modo que la intersección será la media de este grupo de referencia. Los demás coeficientes serán las diferencias entre los otros grupos y el grupo de referencia.
cat_lm <-lm(Y ~ treatment)summary(cat_lm)
Call:
lm(formula = Y ~ treatment)
Residuals:
Min 1Q Median 3Q Max
-5.4227 -1.2765 0.1768 1.3528 4.2373
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.8843 0.2935 9.827 < 2e-16 ***
treatmentlow -1.7197 0.4151 -4.143 6.51e-05 ***
treatmentmed -0.9119 0.4151 -2.197 0.03 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.856 on 117 degrees of freedom
Multiple R-squared: 0.1281, Adjusted R-squared: 0.1132
F-statistic: 8.592 on 2 and 117 DF, p-value: 0.0003298
Entonces, el grupo “high” tiene una media de 2.65, y la diferencia entre las medias del grupo “low” y el grupo “high” es de -0.66, y la diferencia entre el grupo “med” y el grupo “high” es de -1.48. Si deseas tener otro grupo como grupo de referencia, puedes usar relevel para recodificar tu factor de tratamiento.
Call:
lm(formula = Y ~ treatment)
Residuals:
Min 1Q Median 3Q Max
-5.4227 -1.2765 0.1768 1.3528 4.2373
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.1647 0.2935 3.968 0.000125 ***
treatmenthigh 1.7197 0.4151 4.143 6.51e-05 ***
treatmentmed 0.8077 0.4151 1.946 0.054057 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.856 on 117 degrees of freedom
Multiple R-squared: 0.1281, Adjusted R-squared: 0.1132
F-statistic: 8.592 on 2 and 117 DF, p-value: 0.0003298
boxplot(Y ~ treatment)
Ahora la intersección es la media del grupo “low”, y todos los demás coeficientes son las diferencias entre el grupo “low” y los demás. Otra cosa que puedes hacer es poner -1 en el modelo para eliminar la intersección y tener solo las medias de cada grupo como coeficientes.
cat_lm <-lm(Y ~ treatment -1)summary(cat_lm)
Call:
lm(formula = Y ~ treatment - 1)
Residuals:
Min 1Q Median 3Q Max
-5.4227 -1.2765 0.1768 1.3528 4.2373
Coefficients:
Estimate Std. Error t value Pr(>|t|)
treatmentlow 1.1647 0.2935 3.968 0.000125 ***
treatmenthigh 2.8843 0.2935 9.827 < 2e-16 ***
treatmentmed 1.9724 0.2935 6.720 6.96e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.856 on 117 degrees of freedom
Multiple R-squared: 0.5737, Adjusted R-squared: 0.5628
F-statistic: 52.49 on 3 and 117 DF, p-value: < 2.2e-16
Ahora, los tres coeficientes son las medias de los grupos.
Contraste de los coeficientes También podemos ver los contrastes; estos son las diferencias entre todos los pares de grupos. Carga el paquete multcomp y utiliza glht (hipótesis lineales generales) para examinar todas las diferencias de pares.
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = Y ~ treatment - 1)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
high - low == 0 1.7197 0.4151 4.143 0.000186 ***
med - low == 0 0.8077 0.4151 1.946 0.130504
med - high == 0 -0.9119 0.4151 -2.197 0.075894 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
Cada línea de esta salida compara dos grupos entre sí. La primera línea, por ejemplo, compara el grupo “high” con el grupo “low”. Por lo tanto, la diferencia entre las medias de los grupos “high” y “low” es de 1.84. Los valores de p y los intervalos de confianza proporcionados por glht controlan las pruebas múltiples, lo cual es útil. Si deseas ver los intervalos de confianza para las diferencias entre los grupos.
confint(cont)
Simultaneous Confidence Intervals
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = Y ~ treatment - 1)
Quantile = 2.3741
95% family-wise confidence level
Linear Hypotheses:
Estimate lwr upr
high - low == 0 1.71966 0.73419 2.70512
med - low == 0 0.80774 -0.17773 1.79320
med - high == 0 -0.91192 -1.89739 0.07354
Nota: En un modelo con múltiples covariables, las mismas reglas siguen aplicándose en cuanto a las interpretaciones condicionales y marginales de los coeficientes.
Interpretación de los coeficientes en modelos lineales generalizados En los modelos lineales, la interpretación de los parámetros del modelo es lineal, como se discutió anteriormente. Para los modelos lineales generalizados, ahora lee la página de tutoriales sobre la interpretación de los coeficientes en esos modelos.